#region License
// Copyright 2007-2008 Rory Plaire (codekaizen@gmail.com)
/*
 * This license governs use of the accompanying software. If you use the software, you
 * accept this license. If you do not accept the license, do not use the software.
 * 
 * 1. Definitions
 * The terms "reproduce," "reproduction," "derivative works," and "distribution" have the
 * same meaning here as under U.S. copyright law.
 * A "contribution" is the original software, or any additions or changes to the software.
 * A "contributor" is any person that distributes its contribution under this license.
 * "Licensed patents" are a contributor's patent claims that read directly on its contribution.
 * 
 * 2. Grant of Rights
 * (A) Copyright Grant- Subject to the terms of this license, including the license conditions 
 *     and limitations in section 3, each contributor grants you a non-exclusive, worldwide, 
 *     royalty-free copyright license to reproduce its contribution, prepare derivative works 
 *     of its contribution, and distribute its contribution or any derivative works that you create.
 * (B) Patent Grant- Subject to the terms of this license, including the license conditions and 
 *     limitations in section 3, each contributor grants you a non-exclusive, worldwide, royalty-free 
 *     license under its licensed patents to make, have made, use, sell, offer for sale, import, 
 *     and/or otherwise dispose of its contribution in the software or derivative works of the 
 *     contribution in the software.
 * 
 * 3. Conditions and Limitations
 * (A) No Trademark License- This license does not grant you rights to use any contributors' 
 *     name, logo, or trademarks.
 * (B) If you bring a patent claim against any contributor over patents that you claim are 
 *     infringed by the software, your patent license from such contributor to the software 
 *     ends automatically.
 * (C) If you distribute any portion of the software, you must retain all copyright, patent, 
 *     trademark, and attribution notices that are present in the software.
 * (D) If you distribute any portion of the software in source code form, you may do so only under 
 *     this license by including a complete copy of this license with your distribution. If you 
 *     distribute any portion of the software in compiled or object code form, you may only do so 
 *     under a license that complies with this license.
 * (E) The software is licensed "as-is." You bear the risk of using it. The contributors give no 
 *     express warranties, guarantees or conditions. You may have additional consumer rights under 
 *     your local laws which this license cannot change. To the extent permitted under your local laws, 
 *     the contributors exclude the implied warranties of merchantability, fitness for a particular 
 *     purpose and non-infringement.
 * 
 */
#endregion

// Based initially on JAMA : A Java Matrix Package by 
// National Institute of Standards and Technology (NIST) 
// and The MathWorks, a public-domain work.

using System;
using System.Collections.Generic;
using System.Text;
using NPack.Interfaces;

namespace NPack.Matrix
{
    /// <summary>
    /// A matrix. The fundamental unit of computation in linear algebra.
    /// </summary>
    /// <typeparam name="T">Type of element in the matrix.</typeparam>
    public class Matrix<T> : ITransformMatrix<T>, ITransformMatrix<T, Vector<T>, Matrix<T>>
        where T : IEquatable<T>, IComparable<T>, IComputable<Double, T>, IConvertible, IFormattable
    {
        private static readonly MersenneTwister _randomGenerator = new MersenneTwister();

        internal readonly T[] _elements;
        private readonly Int32 _rowCount;
        private readonly Int32 _columnCount;
        private readonly MatrixFormat _format;

        #region Matrix Constructors

        /// <summary>
        /// Creates a zero rectangular Matrix{T} of the given <paramref name="rowCount"/> by <paramref name="columnCount"/>.
        /// </summary>
        /// <param name="format">The format of the matrix, either row-major or column-major.</param>
        /// <param name="rowCount">Number of rows.</param>
        /// <param name="columnCount">Number of columns.</param>
        /// <exception cref="ArgumentOutOfRangeException">Thrown when rowCount or columnCount is less than 1.</exception>
        public Matrix(MatrixFormat format, Int32 rowCount, Int32 columnCount)
        {
            if (rowCount <= 0)
            {
                throw new ArgumentOutOfRangeException("rowCount", rowCount, "Cannot create a matrix without rows.");
            }

            if (columnCount <= 0)
            {
                throw new ArgumentOutOfRangeException("columnCount", columnCount,
                                                      "Cannot create a matrix without columns.");
            }

            _rowCount = rowCount;
            _columnCount = columnCount;

            _format = format;
            _elements = new T[rowCount * columnCount];
        }

        /// <summary>
        /// Creates a rectangular Matrix{T} of the given <paramref name="rowCount"/> 
        /// by <paramref name="columnCount"/> with 
        /// <paramref name="value"/> assigned to the diagonal.
        /// </summary>
        /// <param name="format">The format of the matrix, either row-major or column-major.</param>
        /// <param name="rowCount">Number of rows.</param>
        /// <param name="columnCount">Number of columns.</param>
        /// <param name="value">The value to assign to the diagonal.</param>
        /// <exception cref="ArgumentOutOfRangeException">Thrown when rowCount or columnCount is less than 1.</exception>
        public Matrix(MatrixFormat format, Int32 rowCount, Int32 columnCount, T value)
            : this(format, rowCount, columnCount)
        {
            for (Int32 i = 0; i < rowCount; i++)
            {
                _elements[i * rowCount + i] = value;
            }
        }

        /// <summary>
        /// Creates a new matrix with the given <paramref name="elements"/>.
        /// </summary>
        /// <param name="format">The format of the matrix, either row-major or column-major.</param>
        /// <param name="elements">The elements to use to fill the matrix.</param>
        /// <exception cref="ArgumentNullException">Thrown when parameter 'elements' is null.</exception>
        /// <exception cref="ArgumentOutOfRangeException">Thrown when length of 'elements' is 0 
        /// or when the first element of 'elements' is an array of length 0.</exception>
        /// <exception cref="ArgumentException">When the arrays of 'elements' are not all the same length.</exception>
        public Matrix(MatrixFormat format, T[][] elements)
        {
            if (elements == null)
            {
                throw new ArgumentNullException("elements");
            }

            if (elements.Length == 0 || elements[0] == null || elements[0].Length == 0)
            {
                throw new ArgumentException("Must have at least one element.", "elements");
            }

            _format = format;

            _rowCount = format == MatrixFormat.RowMajor ? elements.Length : elements[0].Length;
            _columnCount = format == MatrixFormat.ColumnMajor ? elements.Length : elements[0].Length;

            _elements = new T[_rowCount * _columnCount];

            for (Int32 i = 0; i < (format == MatrixFormat.RowMajor ? _rowCount : _columnCount); i++)
            {
                if (_columnCount != elements[i].Length)
                {
                    throw new ArgumentException(
                        "Cannot create a matrix which has rows having different numbers of columns.");
                }

                for (Int32 j = 0; j < (format == MatrixFormat.ColumnMajor ? _rowCount : _columnCount); j++)
                {
                    this[i, j] = elements[i][j];
                }
            }
        }

        /// <summary>
        /// Creates a new matrix as a copy of the given matrix.
        /// </summary>
        /// <param name="matrixToCopy">The matrix to copy.</param>
        public Matrix(Matrix<T> matrixToCopy)
            : this(matrixToCopy.Format, matrixToCopy.RowCount, matrixToCopy.ColumnCount)
        {
            MatrixProcessor<T>.SetMatrix(matrixToCopy, this);
        }

        #endregion

        #region GetHashCode

        public override Int32 GetHashCode()
        {
            return (Int32)Math.Pow(RowCount, ColumnCount);
        }

        #endregion

        #region ToString

        public override String ToString()
        {
            StringBuilder buffer = new StringBuilder();

            for (Int32 i = 0; i < RowCount; i++)
            {
                buffer.Append("\n[ ");

                for (Int32 j = 0; j < ColumnCount; j++)
                {
                    buffer.AppendFormat(System.Globalization.CultureInfo.InvariantCulture, "{0:N}, ", this[i, j]);
                }

                buffer.Length -= 2;
                buffer.Append(" ]  ");
            }

            buffer.Length -= 2;

            return buffer.ToString();
        }

        #endregion

        #region Equality Computation

        public override Boolean Equals(Object obj)
        {
            return Equals(obj as Matrix<T>);
        }

        public Boolean Equals(IMatrix<T> other)
        {
            if (ReferenceEquals(this, other))
            {
                return true;
            }

            if (other == null)
            {
                return false;
            }

            if (RowCount != other.RowCount || ColumnCount != other.ColumnCount)
            {
                return false;
            }

            for (Int32 i = 0; i < RowCount; i++)
            {
                for (Int32 j = 0; j < ColumnCount; j++)
                {
                    if (!this[i, j].Equals(other[i, j]))
                    {
                        return false;
                    }
                }
            }

            return true;
        }

        /// <summary>
        /// Compares two <see cref="Matrix{T}"/> instances for element-by-element
        /// equality.
        /// </summary>
        /// <param name="left">The left hand side of the equality comparison.</param>
        /// <param name="right">The right hand side of the equality comparison.</param>
        /// <returns>
        /// <see langword="true"/> if the matrixes are element-by-element equal, <see langword="false"/> otherwise.
        /// </returns>
        public static Boolean operator ==(Matrix<T> left, Matrix<T> right)
        {
            return Equals(left, right);
        }

        /// <summary>
        /// Compares two <see cref="Matrix{T}"/> instances element-by-element for 
        /// inequality.
        /// </summary>
        /// <param name="left">The left hand side of the inequality comparison.</param>
        /// <param name="right">The right hand side of the inequality comparison.</param>
        /// <returns>
        /// <see langword="true"/> if the matrixes are not equal in an element-by-element comparison, 
        /// <see langword="false"/> if they are equal.
        /// </returns>
        public static Boolean operator !=(Matrix<T> left, Matrix<T> right)
        {
            return !Equals(left, right);
        }

        #endregion

        #region Element Access

        /// <summary>
        /// Indexer for the matrix.
        /// </summary>
        /// <param name="row">The row to access.</param>
        /// <param name="column">The column to access.</param>
        /// <returns>The value at the given (row, column) of the matrix.</returns>
        public T this[Int32 row, Int32 column]
        {
            get
            {
                if (row < 0 || row >= RowCount)
                {
                    throw new ArgumentOutOfRangeException("row", row, "Row indexer must be between 0 and RowCount.");
                }

                if (column < 0 || column >= ColumnCount)
                {
                    throw new ArgumentOutOfRangeException("column", column,
                                                          "Column indexer must be between 0 and ColumnCount.");
                }

                Int32 index = computeIndex(row, column);

                return _elements[index];
            }
            set
            {
                if (row < 0 || row >= RowCount)
                {
                    throw new ArgumentOutOfRangeException("row", row, "Row indexer must be between 0 and RowCount.");
                }

                if (column < 0 || column >= ColumnCount)
                {
                    throw new ArgumentOutOfRangeException("column", column,
                                                          "Column indexer must be between 0 and ColumnCount.");
                }

                Int32 index = computeIndex(row, column);

                _elements[index] = value;
            }
        }

        #endregion

        #region Identifying Matrix Properties

        /// <summary>
        /// Gets the number of rows in the matrix.
        /// </summary>
        public Int32 RowCount
        {
            get { return _rowCount; }
        }

        /// <summary>
        /// Gets the number of columns in the matrix.
        /// </summary>
        public Int32 ColumnCount
        {
            get { return _columnCount; }
        }

        /// <summary>
        /// Gets <see langword="true"/> if the <see cref="Matrix{T}"/> is square, false if not.
        /// </summary>
        public Boolean IsSquare
        {
            get { return (RowCount == ColumnCount); }
        }

        /// <summary>
        /// Gets <see langword="true"/> if the <see cref="Matrix{T}"/> is symmetrical, false if not.
        /// </summary>
        public Boolean IsSymmetrical
        {
            get
            {
                if (IsSquare)
                {
                    for (Int32 i = 0; i < _rowCount; i++)
                    {
                        for (Int32 j = 0; j <= i; j++)
                        {
                            if (!this[i, j].Equals(this[j, i]))
                            {
                                return false;
                            }
                        }
                    }

                    return true;
                }

                return false;
            }
        }

        /// <summary>
        /// Gets true if the matrix is singular (non-invertible).
        /// </summary>
        public Boolean IsSingular
        {
            get { return !MatrixProcessor<T>.GetLuDecomposition(this).IsNonsingular; }
        }

        /// <summary>
        /// Gets true if the matrix is invertible (non-singular).
        /// </summary>
        public Boolean IsInvertible
        {
            get { return !IsSingular; }
        }

        /// <summary>
        /// Gets the determinant if the <see cref="Matrix{T}"/> instance is square.
        /// </summary>
        /// <exception cref="InvalidOperationException">Thrown if matrix is not sqare.</exception>
        public Double Determinant
        {
            get { return MatrixProcessor<T>.Determinant(this); }
        }

        /// <summary>
        /// Gets a <see cref="MatrixFormat"/> value indicating how the elements are stored.
        /// </summary>
        public MatrixFormat Format
        {
            get { return _format; }
        }

        #endregion

        #region Matrix Operations

        /// <summary>
        /// Gets a submatrix.
        /// </summary>
        /// <param name="i0">Initial row index.</param>
        /// <param name="i1">Final row index.</param>
        /// <param name="j0">Initial column index.</param>
        /// <param name="j1">Final column index.</param>
        /// <returns>
        /// A submatrix with rows given by the rows <paramref name="i0" />
        /// through <paramref name="i1"/>
        /// and columns <paramref name="j0"/> through <paramref name="j1"/>.
        /// </returns>
        public Matrix<T> GetMatrix(Int32 i0, Int32 i1, Int32 j0, Int32 j1)
        {
            Matrix<T> x = new Matrix<T>(Format, i1 - i0 + 1, j1 - j0 + 1);

            for (Int32 i = i0; i <= i1; i++)
            {
                for (Int32 j = j0; j <= j1; j++)
                {
                    x[i - i0, j - j0] = this[i, j];
                }
            }

            return x;
        }

        /// <summary>
        /// Gets a submatrix.
        /// </summary>
        /// <param name="rowIndicies">The indexes of the rows to include.</param>
        /// <param name="j0">The starting column to include.</param>
        /// <param name="j1">The ending column to include.</param>
        /// <returns>
        /// A submatrix with rows given by <paramref name="rowIndicies"/> 
        /// and columns <paramref name="j0"/> through <paramref name="j1"/>.
        /// </returns>
        public Matrix<T> GetMatrix(Int32[] rowIndicies, Int32 j0, Int32 j1)
        {
            if (j0 > j1)
            {
                throw new ArgumentException("Column index j0 must be less than column index j1");
            }

            if (j0 < 0 || j0 >= ColumnCount)
            {
                throw new ArgumentOutOfRangeException("j0", j0,
                                                      "Parameter must be at least 0 and less than matrix column count.");
            }

            if (j1 < 0 || j1 >= ColumnCount)
            {
                throw new ArgumentOutOfRangeException("j1", j1,
                                                      "Parameter must be at least j0 + 1 and less than matrix column count.");
            }

            Matrix<T> sub = new Matrix<T>(Format, rowIndicies.Length, j1 - j0 + 1);

            for (Int32 i = 0; i < rowIndicies.Length; i++)
            {
                for (Int32 j = j0; j <= j1; j++)
                {
                    if (rowIndicies[i] < 0 || rowIndicies[i] >= RowCount)
                    {
                        throw new ArgumentOutOfRangeException(
                            String.Format("rowIndicies[{0}]", i), rowIndicies[i],
                            "Row index must be at least 0 and less than matrix row count");
                    }

                    sub[i, j - j0] = this[rowIndicies[i], j];
                }
            }

            return sub;
        }

        /// <summary>
        /// Returns a transpose of the <see cref="Matrix{T}"/>.
        /// </summary>
        /// <value>The rows-for-columns, columns-for-rows transpose of A.</value>
        /// <returns>AT, the transpose of matrix A.</returns>
        public Matrix<T> Transpose()
        {
            return MatrixProcessor<T>.Transpose(this);
        }

        /// <summary>
        /// Gets the inverse of the <see cref="Matrix{T}"/> if it is square and non-singular, 
        /// the pseudo-inverse if it is non-square, and null if it is singluar (non-invertible).
        /// </summary>
        public Matrix<T> Inverse
        {
            get { return MatrixProcessor<T>.Invert(this); }
        }

        /// <summary>
        /// Scales the elements in the linear transformation by the given <paramref name="amount"/>.
        /// </summary>
        /// <param name="amount">
        /// Amount of scale to apply uniformly in all dimensions.
        /// </param>
        public void Scale(T amount)
        {
            Scale(amount, MatrixOperationOrder.Default);
        }

        /// <summary>
        /// Scales the elements in the linear transformation by the given <paramref name="amount"/>.
        /// </summary>
        /// <param name="amount">
        /// Amount of scale to apply uniformly in all dimensions.
        /// </param>
        /// <param name="order">The order to apply the transform in.</param>
        public virtual void Scale(T amount, MatrixOperationOrder order)
        {
            Vector<T> scaleVector = new Vector<T>(Format == MatrixFormat.ColumnMajor ? RowCount : ColumnCount);

            for (Int32 i = 0; i < scaleVector.ComponentCount; i++)
            {
                scaleVector[i] = amount;
            }

            Scale(scaleVector, MatrixOperationOrder.Default);
        }

        /// <summary>
        /// Scales the elements in the linear transformation by the given <paramref name="scaleVector"/>.
        /// </summary>
        /// <param name="scaleVector">
        /// Amount of scale to apply on a column-by-column basis.
        /// </param>
        public void Scale(Vector<T> scaleVector)
        {
            Scale(scaleVector, MatrixOperationOrder.Default);
        }

        /// <summary>
        /// Scales the elements in the linear transformation by the given <paramref name="scaleVector"/>.
        /// </summary>
        /// <param name="scaleVector">Amount of scale to apply on a column-by-column basis.</param>
        /// <param name="order">The order to apply the transform in.</param>
        public virtual void Scale(Vector<T> scaleVector, MatrixOperationOrder order)
        {
            Matrix<T> scale = new Matrix<T>(Format, RowCount, ColumnCount, default(T).One);
            MatrixProcessor<T>.Scale(scale, scaleVector);

            Matrix<T> result;

            if (order == MatrixOperationOrder.Prepend)
            {
                if (Format == MatrixFormat.RowMajor)
                {
                    result = MatrixProcessor<T>.Multiply(scale, this);
                }
                else
                {
                    result = MatrixProcessor<T>.Multiply(this, scale);
                }
            }
            else
            {
                if (Format == MatrixFormat.RowMajor)
                {
                    result = MatrixProcessor<T>.Multiply(this, scale);
                }
                else
                {
                    result = MatrixProcessor<T>.Multiply(scale, this);
                }
            }

            MatrixProcessor<T>.SetMatrix(result, this);
        }

        /// <summary>
        /// Applies a shear to the transform by appending the shear to the <see cref="AffineMatrix{T}"/>.
        /// </summary>
        /// <param name="shearVector">The vector used to compute the shear.</param>
        public void Shear(Vector<T> shearVector)
        {
            Shear(shearVector, MatrixOperationOrder.Default);
        }

        /// <summary>
        /// Applies a shear to the transform, either before or after this <see cref="AffineMatrix{T}"/>.
        /// </summary>
        /// <param name="shearVector">The vector used to compute the shear.</param>
        /// <param name="order">The order to apply the transform in.</param>
        public virtual void Shear(Vector<T> shearVector, MatrixOperationOrder order)
        {
            Matrix<T> shear = new Matrix<T>(Format, RowCount, ColumnCount, default(T).One);
            MatrixProcessor<T>.Shear(shear, shearVector);

            Matrix<T> result;

            if (order == MatrixOperationOrder.Prepend)
            {
                if (Format == MatrixFormat.RowMajor)
                {
                    result = MatrixProcessor<T>.Multiply(shear, this);
                }
                else
                {
                    result = MatrixProcessor<T>.Multiply(this, shear);
                }
            }
            else
            {
                if (Format == MatrixFormat.RowMajor)
                {
                    result = MatrixProcessor<T>.Multiply(this, shear);
                }
                else
                {
                    result = MatrixProcessor<T>.Multiply(shear, this);
                }
            }

            MatrixProcessor<T>.SetMatrix(result, this);
        }

        /// <summary>
        /// Rotates the affine transform around the given <paramref name="axis"/>.
        /// </summary>
        /// <param name="axis">The axis to rotate around.</param>
        /// <param name="radians">Angle to rotate through.</param>
        public virtual void RotateAlong(Vector<T> axis, Double radians)
        {
            MatrixProcessor<T>.Rotate(this, axis, radians);
        }

        /// <summary>
        /// Rotates the affine transform around the given <paramref name="axis"/>.
        /// </summary>
        /// <param name="axis">The axis to rotate around.</param>
        /// <param name="radians">Angle to rotate through.</param>
        /// <param name="order">The order to apply the transform in.</param>
        public virtual void RotateAlong(Vector<T> axis, Double radians, MatrixOperationOrder order)
        {
            throw new NotImplementedException();
        }

        /// <summary>
        /// Applies this transform to the given <paramref name="input"/> matrix.
        /// </summary>
        /// <param name="input">Matrix to transform.</param>
        /// <returns>The multiplication of this transform matrix with the input matrix, 
        /// with the transform on the left-hand side of the operation.</returns>
        public Matrix<T> TransformMatrix(Matrix<T> input)
        {
            return Multiply(input);
        }

        ///// <summary>
        ///// Applies this transform to the given <paramref name="input"/> matrix.
        ///// </summary>
        ///// <param name="input">Matrix to transform.</param>
        ///// <returns>The multiplication of this transform matrix with the input matrix, 
        ///// with the transform on the left-hand side of the operation.</returns>
        //public void TransformMatrix(T[][] input)
        //{
        //    MatrixProcessor<T>.Multiply(this, input);
        //}

        /// <summary>
        /// Applies this transform to the given <paramref name="input"/> vector.
        /// </summary>
        /// <param name="input">Vector to transform.</param>
        /// <returns>
        /// The multiplication of this transform matrix with the input vector.
        /// </returns>
        public Vector<T> TransformVector(Vector<T> input)
        {
            T[] elements = new T[input.ComponentCount];
            Array.Copy(input.Components, elements, input.ComponentCount);

            if (Format == MatrixFormat.RowMajor)
            {
                MatrixProcessor<T>.Multiply(elements, this);
            }
            else
            {
                MatrixProcessor<T>.Multiply(this, elements);
            }

            return new Vector<T>(elements);
        }

        /// <summary>
        /// Applies this transform to the given <paramref name="input"/> vectors.
        /// </summary>
        /// <param name="input">Set of vectors to transform.</param>
        public IEnumerable<Vector<T>> TransformVectors(IEnumerable<Vector<T>> input)
        {
            if (Format == MatrixFormat.RowMajor)
            {
                return MatrixProcessor<T>.Multiply(input, this);
            }
            else
            {
                return MatrixProcessor<T>.Multiply(this, input);
            }
        }

        /// <summary>
        /// Applies this transform to the given <paramref name="input"/> vector.
        /// </summary>
        /// <param name="input">Vector to transform.</param>
        public void TransformVector(T[] input)
        {
            if (Format == MatrixFormat.RowMajor)
            {
                MatrixProcessor<T>.Multiply(input, this);
            }
            else
            {
                MatrixProcessor<T>.Multiply(this, input);
            }
        }

        /// <summary>
        /// Applies this transform to the given <paramref name="input"/> vectors.
        /// </summary>
        /// <param name="input">Set of vectors to transform.</param>
        public void TransformVectors(IEnumerable<T[]> input)
        {
            foreach (T[] v in input)
            {
                TransformVector(v);
            }
        }

        #region Arithmetic Operations

        /// <summary>
        /// Returns a copy of the matrix with all the elements negated.
        /// </summary>
        /// <returns>For matrix A, -A.</returns>
        public Matrix<T> Negative()
        {
            return MatrixProcessor<T>.Negate(this);
        }

        /// <summary>
        /// Returns a matrix which is the result of the instance plus the <paramref name="value"/>.
        /// </summary>
        /// <param name="value">The matrix to add.</param>
        /// <returns>For matrix A, A + value.</returns>
        public Matrix<T> Add(Matrix<T> value)
        {
            return MatrixProcessor<T>.Add(this, value);
        }

        /// <summary>
        /// Returns a matrix which is the result of the instance minus the <paramref name="value"/>.
        /// </summary>
        /// <param name="value">The matrix to subtract.</param>
        /// <returns>For matrix A, A - value.</returns>
        public Matrix<T> Subtract(Matrix<T> value)
        {
            return MatrixProcessor<T>.Subtract(this, value);
        }

        /// <summary>
        /// Scalar multiplication.
        /// </summary>
        /// <param name="value">Scalar to multiply.</param>
        /// <returns>For matrix A, (value)A.</returns>
        public Matrix<T> Multiply(T value)
        {
            return MatrixProcessor<T>.ScalarMultiply(this, value);
        }

        /// <summary>
        /// Matrix multiplication.
        /// </summary>
        /// <param name="value">Matrix to multiply.</param>
        /// <returns>For matrix A, (value)A.</returns>
        public Matrix<T> Multiply(Matrix<T> value)
        {
            return MatrixProcessor<T>.Multiply(this, value);
        }

        #region Arithmetical Operators

        /// <summary>
        /// Negates the values of a matrix.
        /// </summary>
        /// <param name="value">The matrix to negate.</param>
        /// <returns>A negated matrix instance.</returns>
        public static Matrix<T> operator -(Matrix<T> value)
        {
            if (value == null)
            {
                throw new ArgumentNullException("value");
            }

            return MatrixProcessor<T>.Negate(value);
        }

        /// <summary>
        /// Adds two matrixes element-by-element.
        /// </summary>
        /// <param name="lhs">The left hand side of the addition operation.</param>
        /// <param name="rhs">The right hand side of the addition operation.</param>
        /// <returns>The sum of <paramref name="lhs"/> and <paramref name="rhs"/>.</returns>
        /// <exception cref="InvalidOperationException">
        /// Thrown if the dimensions of <paramref name="lhs"/> and <paramref name="rhs"/>
        /// are not the same.
        /// </exception>
        public static Matrix<T> operator +(Matrix<T> lhs, Matrix<T> rhs)
        {
            if (lhs == null)
            {
                throw new ArgumentNullException("lhs");
            }

            if (rhs == null)
            {
                throw new ArgumentNullException("rhs");
            }

            return MatrixProcessor<T>.Add(lhs, rhs);
        }

        /// <summary>
        /// Subtracts two matrixes element-by-element.
        /// </summary>
        /// <param name="lhs">The left hand side of the subtraction operation.</param>
        /// <param name="rhs">The right hand side of the subtraction operation.</param>
        /// <returns>The difference of <paramref name="lhs"/> and <paramref name="rhs"/>.</returns>
        /// <exception cref="InvalidOperationException">
        /// Thrown if the dimensions of <paramref name="lhs"/> and <paramref name="rhs"/>
        /// are not the same.
        /// </exception>
        public static Matrix<T> operator -(Matrix<T> lhs, Matrix<T> rhs)
        {
            if (lhs == null)
            {
                throw new ArgumentNullException("lhs");
            }

            if (rhs == null)
            {
                throw new ArgumentNullException("rhs");
            }

            return MatrixProcessor<T>.Subtract(lhs, rhs);
        }

        /// <summary>
        /// Scalar multiplication.
        /// </summary>
        public static Matrix<T> operator *(Matrix<T> lhs, T rhs)
        {
            if (lhs == null)
            {
                throw new ArgumentNullException("lhs");
            }

            return MatrixProcessor<T>.ScalarMultiply(lhs, rhs);
        }

        /// <summary>
        /// <see cref="Matrix{T}"/> multiplication.
        /// </summary>
        public static Matrix<T> operator *(Matrix<T> lhs, Matrix<T> rhs)
        {
            if (lhs == null)
            {
                throw new ArgumentNullException("lhs");
            }

            if (rhs == null)
            {
                throw new ArgumentNullException("rhs");
            }

            return MatrixProcessor<T>.Multiply(lhs, rhs);
        }

        #endregion

        #endregion

        #endregion


        #region IComparable<Matrix<T>> Members

        public Int32 CompareTo(Matrix<T> other)
        {
            if (Equals(other))
                return 0;

            if (_elements.Length > other._elements.Length)
                return 1;
            if (_elements.Length < other._elements.Length)
                return -1;

            if (Format != other.Format)
                return Comparer<MatrixFormat>.Default.Compare(Format, other.Format);

            for (int i = 0; i < _elements.Length; i++)
            {
                int val = Comparer<T>.Default.Compare(_elements[i], other._elements[i]);
                if (val != 0)
                    return val;
            }
            return 0;

        }

        #endregion

        #region IComputable<Matrix<T>> Members

        public Matrix<T> Abs()
        {
            throw new NotImplementedException();
        }

        public Matrix<T> Set(Double value)
        {
            throw new NotImplementedException();
        }

        #endregion

        #region IDivisible<Matrix<T>> Members

        public Matrix<T> Divide(Matrix<T> b)
        {
            throw new NotSupportedException();
        }

        #endregion

        #region IBooleanComparable<Matrix<T>> Members

        public Boolean GreaterThan(Matrix<T> value)
        {
            return CompareTo(value) > 0;
        }

        public Boolean GreaterThanOrEqualTo(Matrix<T> value)
        {
            return CompareTo(value) >= 0;
        }

        public Boolean LessThan(Matrix<T> value)
        {
            return CompareTo(value) < 0;
        }

        public Boolean LessThanOrEqualTo(Matrix<T> value)
        {
            return CompareTo(value) <= 0;
        }

        #endregion

        #region IExponential<Matrix<T>> Members

        public Matrix<T> Power(Double exponent)
        {
            throw new NotImplementedException();
        }

        public Matrix<T> Sqrt()
        {
            throw new NotImplementedException();
        }

        public Matrix<T> Log(Double newBase)
        {
            throw new NotImplementedException();
        }

        public Matrix<T> Log()
        {
            throw new NotImplementedException();
        }

        public Matrix<T> Exp()
        {
            throw new NotImplementedException();
        }

        #endregion

        #region IEquatable<Matrix<T>> Members

        public Boolean Equals(Matrix<T> other)
        {
            if (other == null)
            {
                return false;
            }

            if (other.RowCount != RowCount || other.ColumnCount != ColumnCount)
            {
                return false;
            }

            for (Int32 i = 0; i < RowCount; i++)
            {
                for (Int32 j = 0; j < ColumnCount; j++)
                {
                    if (!other[i, j].Equals(this[i, j]))
                    {
                        return false;
                    }
                }
            }

            return true;
        }

        #endregion

        #region IHasZero<IMatrix<T>> Members

        /// <summary>
        /// Gets a zero matrix of the same row and column rank.
        /// </summary>
        /// <remarks>
        /// Uses <c>default(T)</c> as zero.
        /// </remarks>
        public Matrix<T> Zero
        {
            get { return new Matrix<T>(Format, RowCount, ColumnCount, default(T)); }
        }

        #endregion

        #region IDivisible<IMatrix<T>> Members

        /// <summary>
        /// Matrix division is not defined. Throws a <see cref="NotSupportedException"/>.
        /// </summary>
        /// <exception cref="NotSupportedException"/>
        //IMatrix<T> IDivisible<Matrix<T>>.Divide(IMatrix<T> a)
        //{
        //    throw new NotSupportedException();
        //}

        #endregion

        #region IHasOne<IMatrix<T>> Members

        /// <summary>
        /// Gets an identity matrix with the same row rank.
        /// </summary>
        public Matrix<T> One
        {
            get { return Identity(Format, RowCount); }
        }

        #endregion

        /// <summary>
        /// Creates an element-by-element copy of the matrix.
        /// </summary>
        public Matrix<T> Clone()
        {
            Matrix<T> clone = new Matrix<T>(Format, RowCount, ColumnCount);

            MatrixProcessor<T>.SetMatrix(this, clone);

            return clone;
        }

        #region Static Matrix Generators

        /// <summary>
        /// Generates a square matrix of the given <paramref name="rank"/> with
        /// the number 1 in each element of the diagonal.
        /// </summary>
        /// <param name="format">The format of the matrix, either row-major or column-major.</param>
        /// <param name="rank">Number of rows and columns of the <see cref="Matrix{T}"/>.</param>
        /// <returns>An identiy matrix of the given rank.</returns>
        public static Matrix<T> Identity(MatrixFormat format, Int32 rank)
        {
            return new Matrix<T>(format, rank, rank, default(T).One);
        }

        /// <summary>
        /// Creates a <see cref="Matrix{T}"/> filled with pseudo-random values.
        /// </summary>
        /// <remarks>
        /// The random values are generated by using the Mersenne Twister algorithm.
        /// </remarks>
        /// <param name="format">The format of the matrix, either row-major or column-major.</param>
        /// <param name="columns">Number of columns in the matrix.</param>
        /// <param name="rows">Number of rows in the matrix.</param>
        /// <returns>A <paramref name="rows"/> by <paramref name="columns"/> matrix with random values.</returns>
        public static Matrix<T> CreateRandom(MatrixFormat format, Int32 rows, Int32 columns)
        {
            Matrix<T> random = new Matrix<T>(format, rows, columns);

            for (Int32 i = 0; i < rows; i++)
            {
                for (Int32 j = 0; j < columns; j++)
                {
                    random[i, j] = default(T).Set(_randomGenerator.NextDouble());
                }
            }

            return random;
        }

        #endregion

        private Int32 computeIndex(Int32 row, Int32 column)
        {
            Int32 index;
            index = _format == MatrixFormat.RowMajor
                        ? row * _columnCount + column
                        : column * _rowCount + row;
            return index;
        }

        #region ITransformMatrix<T> Member

        void ITransformMatrix<T>.Scale(T amount)
        {
            throw new NotImplementedException();
        }

        void ITransformMatrix<T>.Scale(T amount, MatrixOperationOrder order)
        {
            throw new NotImplementedException();
        }

        void ITransformMatrix<T>.Scale(IVector<T> scaleVector)
        {
            throw new NotImplementedException();
        }

        void ITransformMatrix<T>.Scale(IVector<T> scaleVector, MatrixOperationOrder order)
        {
            throw new NotImplementedException();
        }

        void ITransformMatrix<T>.Shear(IVector<T> shearVector)
        {
            throw new NotImplementedException();
        }

        void ITransformMatrix<T>.Shear(IVector<T> shearVector, MatrixOperationOrder order)
        {
            throw new NotImplementedException();
        }

        void ITransformMatrix<T>.RotateAlong(IVector<T> axis, double radians)
        {
            throw new NotImplementedException();
        }

        void ITransformMatrix<T>.RotateAlong(IVector<T> axis, double radians, MatrixOperationOrder order)
        {
            throw new NotImplementedException();
        }

        IMatrix<T> ITransformMatrix<T>.TransformMatrix(IMatrix<T> input)
        {
            throw new NotImplementedException();
        }

        IVector<T> ITransformMatrix<T>.TransformVector(IVector<T> input)
        {
            throw new NotImplementedException();
        }

        void ITransformMatrix<T>.TransformVector(T[] input)
        {
            throw new NotImplementedException();
        }

        IEnumerable<IVector<T>> ITransformMatrix<T>.TransformVectors(IEnumerable<IVector<T>> input)
        {
            throw new NotImplementedException();
        }

        void ITransformMatrix<T>.TransformVectors(IEnumerable<T[]> input)
        {
            throw new NotImplementedException();
        }

        #endregion

        #region IMatrix<T> Member

        double IMatrix<T>.Determinant
        {
            get { throw new NotImplementedException(); }
        }

        int IMatrix<T>.ColumnCount
        {
            get { throw new NotImplementedException(); }
        }

        MatrixFormat IMatrix<T>.Format
        {
            get { throw new NotImplementedException(); }
        }

        bool IMatrix<T>.IsSingular
        {
            get { throw new NotImplementedException(); }
        }

        bool IMatrix<T>.IsInvertible
        {
            get { throw new NotImplementedException(); }
        }

        IMatrix<T> IMatrix<T>.Inverse
        {
            get { throw new NotImplementedException(); }
        }

        bool IMatrix<T>.IsSquare
        {
            get { return this.IsSquare; }
        }

        bool IMatrix<T>.IsSymmetrical
        {
            get { return this.IsSymmetrical; }
        }

        int IMatrix<T>.RowCount
        {
            get { return this.RowCount; }
        }

        T IMatrix<T>.this[int row, int column]
        {
            get
            {
                return this[row, column];
            }
            set
            {
                this[row, column] = value;
            }
        }

        IMatrix<T> IMatrix<T>.Clone()
        {
            return this.Clone();
        }

        IMatrix<T> IMatrix<T>.GetMatrix(int[] rowIndexes, int startColumn, int endColumn)
        {
            return this.GetMatrix(rowIndexes, startColumn, endColumn);
        }

        IMatrix<T> IMatrix<T>.Transpose()
        {
            return this.Transpose();
        }

        #endregion

        #region IEquatable<IMatrix<T>> Member

        bool IEquatable<IMatrix<T>>.Equals(IMatrix<T> other)
        {
            return this.Equals((Matrix<T>)other);
        }

        #endregion

        #region IComparable<IMatrix<T>> Member

        int IComparable<IMatrix<T>>.CompareTo(IMatrix<T> other)
        {
            return this.CompareTo((Matrix<T>)other);
        }

        #endregion

        #region IComputable<IMatrix<T>> Member

        IMatrix<T> IComputable<IMatrix<T>>.Abs()
        {
            return this.Abs();
        }

        IMatrix<T> IComputable<IMatrix<T>>.Set(double value)
        {
            return this.Set(value);
        }

        #endregion

        #region INegatable<IMatrix<T>> Member

        IMatrix<T> INegatable<IMatrix<T>>.Negative()
        {
            return this.Negative();
        }

        #endregion

        #region ISubtractable<IMatrix<T>> Member

        IMatrix<T> ISubtractable<IMatrix<T>>.Subtract(IMatrix<T> b)
        {
            return this.Subtract((Matrix<T>)b);
        }

        #endregion

        #region IHasZero<IMatrix<T>> Member

        IMatrix<T> IHasZero<IMatrix<T>>.Zero
        {
            get { throw new NotImplementedException(); }
        }

        #endregion

        #region IAddable<IMatrix<T>> Member

        IMatrix<T> IAddable<IMatrix<T>>.Add(IMatrix<T> b)
        {
            return this.Add((Matrix<T>)b);
        }

        #endregion

        #region IDivisible<IMatrix<T>> Member

        IMatrix<T> IDivisible<IMatrix<T>>.Divide(IMatrix<T> b)
        {
            return this.Divide((Matrix<T>)b);
        }

        #endregion

        #region IHasOne<IMatrix<T>> Member

        IMatrix<T> IHasOne<IMatrix<T>>.One
        {
            get { throw new NotImplementedException(); }
        }

        #endregion

        #region IMultipliable<IMatrix<T>> Member

        IMatrix<T> IMultipliable<IMatrix<T>>.Multiply(IMatrix<T> b)
        {
            return this.Multiply((Matrix<T>)b);
        }

        #endregion

        #region IBooleanComparable<IMatrix<T>> Member

        bool IBooleanComparable<IMatrix<T>>.GreaterThan(IMatrix<T> value)
        {
            return this.GreaterThan((Matrix<T>)value);
        }

        bool IBooleanComparable<IMatrix<T>>.GreaterThanOrEqualTo(IMatrix<T> value)
        {
            return this.GreaterThanOrEqualTo((Matrix<T>)value);
        }

        bool IBooleanComparable<IMatrix<T>>.LessThan(IMatrix<T> value)
        {
            return this.LessThan((Matrix<T>)value);
        }

        bool IBooleanComparable<IMatrix<T>>.LessThanOrEqualTo(IMatrix<T> value)
        {
            return this.LessThanOrEqualTo((Matrix<T>)value);
        }

        #endregion

        #region IExponential<IMatrix<T>> Member

        IMatrix<T> IExponential<IMatrix<T>>.Power(double exponent)
        {
            return this.Power(exponent);
        }

        IMatrix<T> IExponential<IMatrix<T>>.Sqrt()
        {
            return this.Sqrt();
        }

        IMatrix<T> IExponential<IMatrix<T>>.Log(double newBase)
        {
            return this.Log(newBase);
        }

        IMatrix<T> IExponential<IMatrix<T>>.Log()
        {
            return this.Log();
        }

        IMatrix<T> IExponential<IMatrix<T>>.Exp()
        {
            return this.Exp();
        }

        #endregion
    }
}